In [ ]:
using LinearAlgebra
using Symbolics
using Plots
using LaTeXStrings
using Distributions
using MAT
Modelo SIR¶
Modelo SIR¶
In [ ]:
N=1
x,y,z,v_x,v_y,v_z,r0,rs,vs,m,t=symN(N)
γ=1/2
β=1/2
conds=[0 0 0 997 3 0]
NN=sum(conds)
F=[-β*v_x[1]*v_y[1]/NN β*v_x[1]*v_y[1]/NN-γ*v_y[1] γ*v_y[1]]
tf=100
dt=0.01
xs,ys,zs,vxs,vys,vzs,ts=RK5(conds,dt,tf,F,false,vs);
plot(ts,vxs',label="Susceptoibles " *L"S_0="*string(vxs[1]),color="Blue")
plot!(ts,vys',label="Infectados " *L"I_0="*string(vys[1]),color="Red")
plot!(ts,vzs',label="recuperados " *L"R_0="*string(vzs[1]),color="green")
title!("Solucion Modelo SIR "*L" \beta= "*string(β)*" "*L"\gamma="*string(γ) )
Out[ ]:
In [ ]:
γ=1/20
β=1/2
conds=[0 0 0 990 10 0]
NN=sum(conds)
F=[-β*v_x[1]*v_y[1]/NN β*v_x[1]*v_y[1]/NN-γ*v_y[1] γ*v_y[1]]
tf=100
dt=0.01
xs,ys,zs,vxs,vys,vzs,ts=RK5(conds,dt,tf,F,false,vs);
plot(ts,vxs',label="Susceptoibles " *L"S_0="*string(vxs[1]),color="Blue")
plot!(ts,vys',label="Infectados " *L"I_0="*string(vys[1]),color="Red")
plot!(ts,vzs',label="recuperados " *L"R_0="*string(vzs[1]),color="green")
title!("Solucion Modelo SIR "*L" \beta= "*string(β)*" "*L"\gamma="*string(γ) )
Out[ ]:
In [ ]:
β=0.345
γ=0.0743#0.157
conds=[0 0 0 5177 23 0]
NN=sum(conds)
F=[-β*v_x[1]*v_y[1]/NN β*v_x[1]*v_y[1]/NN-γ*v_y[1] γ*v_y[1]]
tf=100
dt=0.01
xs,ys,zs,vxs,vys,vzs,ts=RK5(conds,dt,tf,F,false,vs);
plot(ts,vxs',label="Susceptoibles " *L"S_0="*string(vxs[1]),color="Blue",size=(1000,800))
plot!(ts,vys',label="Infectados " *L"I_0="*string(vys[1]),color="Red")
plot!(ts,vzs',label="recuperados " *L"R_0="*string(vzs[1]),color="green")
title!("Solucion Modelo SIR "*L" \beta= "*string(β)*" "*L"\gamma="*string(γ) )
Out[ ]:
In [ ]:
Sr,Ir,Rr=ruido(vxs',vys',vzs',0.0);
# Create a sample matrix
#data = Dict("S" => Sr, "I" => Ir, "R" => Rr, "t" => ts,"N"=>NN)
data = Dict("S" => Sr[1,1:100:end], "I" => Ir[1:100:end], "R" => Rr[1:100:end], "t" => ts[1:100:end],"N"=>NN)
# Specify the file name (change this to your desired file name)
filename = "SIR_limpio.mat"
# Save the data to the .mat file
MAT.matwrite(filename, data)
In [ ]:
Sr,Ir,Rr=ruido(vxs',vys',vzs',0.05);
# Create a sample matrix
#data = Dict("S" => Sr, "I" => Ir, "R" => Rr, "t" => ts,"N"=>NN)
data = Dict("S" => Sr[1:100:end], "I" => Ir[1:100:end], "R" => Rr[1:100:end], "t" => ts[1:100:end],"N"=>NN)
# Specify the file name (change this to your desired file name)
filename = "SIR_ejemp2.mat"
# Save the data to the .mat file
MAT.matwrite(filename, data)
In [ ]:
Is[1]
22.86423174304127
In [ ]:
data = matread("C:\\Users\\Arif\\OneDrive - Instituto Tecnologico y de Estudios Superiores de Monterrey\\Estocasticos\\Reto\\Fase 3\\SIRs1.mat")
Ss=data["Ss"];
Is=data["Is"];
Rs=data["Rs"];
tss=data["ts"];
plot!(tss,Ss',label="Susceptoibles Suavizados " *L"S_0="*string(Ss[1]),color="Blue",line = (:dash, 2))
plot!(tss,Is',label="Infectados Suavizados " *L"I_0="*string(Is[1]),color="Red",line = (:dash, 2))
plot!(tss,Rs',label="recuperados Suavizados" *L"R_0="*string(Rs[1]),color="green",line = (:dash, 2))
Out[ ]:
Pregunta 4 Fase 3¶
In [ ]:
β=0.5785#0.0079
γ=0.0597#0.0008#0.157
conds=[0 0 0 1223 7 0]
NN=sum(conds)
F=[-β*v_x[1]*v_y[1]/NN β*v_x[1]*v_y[1]/NN-γ*v_y[1] γ*v_y[1]]
tf=50
dt=0.01
xs,ys,zs,vxs,vys,vzs,ts=RK5(conds,dt,tf,F,false,vs);
plot(ts,vxs',label="Susceptoibles " *L"S_0="*string(vxs[1]),color="Blue",size=(1000,800))
plot!(ts,vys',label="Infectados " *L"I_0="*string(vys[1]),color="Red")
plot!(ts,vzs',label="recuperados " *L"R_0="*string(vzs[1]),color="green")
title!("Solucion Modelo SIR "*L" \beta= "*string(β)*" "*L"\gamma="*string(γ) )
Out[ ]:
In [ ]:
Is[1]
Out[ ]:
9.040857028704194
In [ ]:
filename = "C:\\Users\\Arif\\OneDrive - Instituto Tecnologico y de Estudios Superiores de Monterrey\\Estocasticos\\Reto\\Fase 3\\matriz_combinada.mat"
mat_data = matread(filename)
Ss = mat_data["matriz_combinada"][1, :]
Is = mat_data["matriz_combinada"][2, :]
Rs = mat_data["matriz_combinada"][3, :]
tss = mat_data["matriz_combinada"][4, :]
plot!(tss,Ss,label="Susceptoibles Suavizados " *L"S_0="*string(Ss[1]),color="Blue",line = (:dash, 2))
plot!(tss,Is,label="Infectados Suavizados " *L"I_0="*string(Is[1]),color="Red",line = (:dash, 2))
plot!(tss,Rs,label="recuperados Suavizados" *L"R_0="*string(Rs[1]),color="green",line = (:dash, 2))
Out[ ]:
In [ ]:
maximum(tss)
100.0
In [ ]:
plot(tss[1:end-1],diff(Ss,dims=2)',label="Suavizado")
plot!(tss[1:end],-(0.345.*(Is.*Ss.*diff(tss,dims=1)[1])./NN)',label="optimizado")
In [ ]:
(Is.*Ss.*diff(tss,dims=1)[1])./N
1×100 Matrix{Float64}:
1.18372e5 1.23605e5 1.3212e5 … 16261.0 15651.8 15071.7 14531.3
In [ ]:
γ=0.1
β=0.2
conds=[0 0 0 998 1 0]
NN=sum(conds)
F=[-β*v_x[1]*v_y[1]/NN β*v_x[1]*v_y[1]/NN-γ*v_y[1] γ*v_y[1]]
tf=100
dt=0.01
xs,ys,zs,vxs,vys,vzs,ts=RK5(conds,dt,tf,F,false,vs);
plot(ts,vxs',label="Susceptoibles " *L"S_0="*string(vxs[1]),color="Blue")
plot!(ts,vys',label="Infectados " *L"I_0="*string(vys[1]),color="Red")
plot!(ts,vzs',label="recuperados " *L"R_0="*string(vzs[1]),color="green")
title!("Solucion Modelo SIR "*L" \beta= "*string(β)*" "*L"\gamma="*string(γ) )
Out[ ]:
In [ ]:
F'
Out[ ]:
$$ \begin{equation}
\left[
\begin{array}{c}
- 0.0002002 v_{x_1} v_{y_1} \\
- 0.1 v_{y_1} + 0.0002002 v_{x_1} v_{y_1} \\
0.1 v_{y_1} \\
\end{array}
\right]
\end{equation}
$$
Plot N -parts¶
In [ ]:
function plotN(xs,ys,zs,N)
gr(size=(800,800))
if length(xs[1,:])>=1000
xs=xs[:,1:100:end]
ys=ys[:,1:100:end]
zs=zs[:,1:100:end]
end
#plot_layout = @layout [a,b]
p = plot()#plot(layout=plot_layout)
if length(size(xs))>1
N=size(xs)[1]
else
N=1
end
if N==1
scatter!([xs[N,end]],[ys[N,end]],[zs[N,end]],label="Posición Final particula "*string(N),color="red",markersize=4,aspect_ratio= :equal,xlabel=L"X \mu m",ylabel=L"Y \mu m")
plot!(xs[:],ys[:],zs[:],label=false,color="red",lw=2,alpha=0.7) #
else
for i in 1:N
#if i>1
#plot(xs[i,:],ys[i,:],label=false,color=i,lw=2,xlim = [minimum(xs), maximum(xs)],ylim = [minimum(ys), maximum(ys)])
#end
plot!(xs[i,:],ys[i,:],zs[i,:],label=false,color=i,lw=2,xlim = [minimum(xs), maximum(xs)],ylim = [minimum(ys), maximum(ys)],xlabel=L"X \mu m",ylabel=L"Y \mu m",aspect_ratio=:equal,alpha=0.7)
scatter!([xs[i,end]],[ys[i,end]],[zs[i,end]],label="Posición Final particula "*string(i),color=i,markersize=4)
end
end
display(p)
end
Out[ ]:
plotN (generic function with 1 method)
Runge Kutta 5to Orden¶
In [ ]:
function RK5(c0,dt,tf,F1,R,vs)
F1=hcat(vs,F1,zeros(N))
F=build_function(F1,[x...]...,[y...]...,[z...]...,[v_x...]...,[v_y...]...,[v_z...]...,t,expression=Val{false})
##vectores
ts=collect(0:dt:tf)
n=length(ts)
####################################Posiciones
xs=zeros(N,n)
ys=zeros(N,n)
zs=zeros(N,n)
vxs=zeros(N,n)
vys=zeros(N,n)
vzs=zeros(N,n)
#####################################Velocidades
#############Condiciones iniciales
xs[:,1]=c0[:,1]
ys[:,1]=c0[:,2]
zs[:,1]=c0[:,3]
vxs[:,1]=c0[:,4]
vys[:,1]=c0[:,5]
vzs[:,1]=c0[:,6]
for i in 1:n-1
#Runge Kutta
z0=hcat(xs[:,i],ys[:,i],zs[:,i],vxs[:,i],vys[:,i],vzs[:,i],ones(N).*ts[i])
k1=F[1](z0...)
k2=F[1](z0.+(dt*k1)/2...)
k3=F[1](z0.+dt*(3*k1+k2)/16...)
k4=F[1](z0.+dt*k3/2...)
k5=F[1](z0.+dt*(-3*k2+6*k3+9*k4)/16...)
k6=F[1](z0.+dt*(k1+4*k2+6*k3-12*k4+8*k5)/7...)
k1=7*(k1+k6)
k3=16*k3
k2=16*k5
k4=12*k4
xs[:,i+1]=xs[:,i]+(dt/90)*(k1[:,1]+2*k2[:,1]+2*k3[:,1]+k4[:,1])
ys[:,i+1]=ys[:,i]+(dt/90)*(k1[:,2]+2*k2[:,2]+2*k3[:,2]+k4[:,2])
zs[:,i+1]=zs[:,i]+(dt/90)*(k1[:,3]+2*k2[:,3]+2*k3[:,3]+k4[:,3])
vxs[:,i+1]=vxs[:,i]+(dt/90)*(k1[:,4]+2*k2[:,4]+2*k3[:,4]+k4[:,4])
vys[:,i+1]=vys[:,i]+(dt/90)*(k1[:,5]+2*k2[:,5]+2*k3[:,5]+k4[:,5])
vzs[:,i+1]=vzs[:,i]+(dt/90)*(k1[:,6]+2*k2[:,6]+2*k3[:,6]+k4[:,6])
#axs3[i]=h1(xs1[i],xs2[i],xs3[i],ys1[i],ys2[i],ys3[i],zs1[i],zs2[i],zs3[i],vxs1[i],vxs2[i],vxs3[i],vys1[i],vys2[i],vys3[i],vzs1[i],vzs2[i],vzs3[i],ts[i])
#ays3[i]=h2(xs1[i],xs2[i],xs3[i],ys1[i],ys2[i],ys3[i],zs1[i],zs2[i],zs3[i],vxs1[i],vxs2[i],vxs3[i],vys1[i],vys2[i],vys3[i],vzs1[i],vzs2[i],vzs3[i],ts[i])
#azs3[i]=h3(xs1[i],xs2[i],xs3[i],ys1[i],ys2[i],ys3[i],zs1[i],zs2[i],zs3[i],vxs1[i],vxs2[i],vxs3[i],vys1[i],vys2[i],vys3[i],vzs1[i],vzs2[i],vzs3[i],ts[i]);
end
if R==true
nn=Int(round(n/1000))
gr(size=(800,500))
plot(xs[1,1:nn:end],ys[1,1:nn:end],zs[1,1:nn:end],label="masa 1",color="yellow",lw=5)
#plot!(xs2[1:nn:end],ys2[1:nn:end],zs2[1:nn:end],label="masa 2",color="blue",lw=2)
#plot!(xs3[1:nn:end],ys3[1:nn:end],zs3[1:nn:end],label="masa 3",color="red",lw=2)
end
return xs,ys,zs,vxs,vys,vzs,ts
########Graficas####################################
end
Out[ ]:
RK5 (generic function with 1 method)
Animacion 3D¶
In [ ]:
#xs1=
n=length(ts)
mm=5
#xx=vcat(xs1,xs2,xs3)
#yy=vcat(ys1,ys2,ys3)
#zz=vcat(zs1,zs2,zs3)
ext=1/3
l=1 #rotacion de la grafica
anim3d=@animate for i in 1:10:n;
gr(size=(1000,800))
a=1
j=1
if i>1
j=Int(ceil(i*0.03))
a=range(0,1,length=(i-j+1))
end
Plots.plot(xs[1,j:i],ys[1,j:i],zs[1,j:i],color="yellow",lw=2
#Plots.plot(xs1[1:i],ys1[1:i],zs1[1:i],color="yellow",lw=2,xlim = [minimum(xx)-ext, maximum(xx)+ext],
,label="",xlabel="x",ylabel="y",zlabel="z",camera = (i/4*l, 20)
,bgcolor="black",grid=false,alpha=a);
Plots.scatter!([xs[1,i]],[ys[1,i]],[zs[1,i]], color = "yellow", markersize = mm,label="Masa1");
#Plots.plot!(xs[2,j:i],ys[2,j:i],zs[2,j:i],color="blue",lw=2,label="",alpha=a);
#Plots.scatter!([xs[2,i]],[ys[2,i]],[zs[2,i]], color = "blue", markersize = mm ,label="Masa2");
#Plots.plot!(xs[3,j:i],ys[3,j:i],zs[3,j:i],color="red",lw=2,label="",alpha=a);
#Plots.scatter!([xs[3,i]],[ys[3,i]],[zs[3,i]], color = "red", markersize = mm,label="Masa3",markershape=:star4);
end;
In [ ]:
gif(anim3d,"jimy .gif",fps=24)
[ Info: Saved animation to C:\Users\Arif\OneDrive - Instituto Tecnologico y de Estudios Superiores de Monterrey\Escuela\Carrera\4to Semestre\mecánica clásica\Reto\jimy .gif
Out[ ]:
Animacion 2D¶
In [ ]:
n=length(ts)
mm=5
xx=vcat(xs1,xs2,xs3)
yy=vcat(ys1,ys2,ys3)
sc=1/100
anim2d=@animate for i in 1:10:n;
gr(size=(1000,1000))
Legend=(bgcolor=(:dark))
a=1
if i>1
a=range(0,1,length=i)
end
Plots.plot(xs1[1:i],ys1[1:i],color="yellow",lw=2,xlim = [minimum(xx)-1, maximum(xx)+1],
ylim = [minimum(yy)-1, maximum(yy)+1],label="",xlabel="x",ylabel="y",bgcolor="black",grid=false,axis=false,alpha=a);
Plots.scatter!([xs1[i]],[ys1[i]], color = "yellow", markersize = mm,label="Masa1");
Plots.quiver!([xs1[i]],[ys1[i]],quiver=([axs1[i].*sc],[ays1[i].*sc]),color="yellow",lw=3,length=1/80)
Plots.plot!(xs2[1:i],ys2[1:i],color="blue",lw=2,label="",alpha=a);
#Plots.scatter!([xs2[i]],[ys2[i]], color = "blue", markersize = mm ,label="Masa2");
#Plots.quiver!([xs2[i]],[ys2[i]],quiver=([axs2[i].*sc],[ays2[i].*sc]),color="blue",lw=3,length=1/80)
Plots.plot!(xs3[1:i],ys3[1:i],color="red",lw=2,label="",alpha=a);
#Plots.scatter!([xs3[i]],[ys3[i]], color = "red", markersize = mm,label="Masa3");
#Plots.quiver!([xs3[i]],[ys3[i]],quiver=([axs3[i].*sc],[ays3[i].*sc]),color="red",lw=3,length=1/80,normalize=true)
#
end;
In [ ]:
gif(anim2d,"Ciclotron2d.gif",fps=24)
[ Info: Saved animation to C:\Users\Arif\OneDrive - Instituto Tecnologico y de Estudios Superiores de Monterrey\Escuela\Carrera\4to Semestre\mecánica clásica\Reto\Ciclotron2d.gif
Out[ ]:
Animación¶
In [ ]:
function animac(xs,ys,zs,N,n)
anim2d=@animate for i in 1:50:n;
gr(size=(1000,1000))
a=1
j=1
if i>1
j=Int(ceil(i*0.3))
a=range(0,1,length=(i-j+1))
end
for d in 1:N
if d==1
Plots.plot(xs[d,j:i],ys[d,j:i],zs[d,j:i],color=d,lw=2,label="",xlabel=L"X \mu m",ylabel=L"Y \mu m",bgcolor="black",alpha=a,xlim = [minimum(xs), maximum(xs)],ylim = [minimum(ys), maximum(ys)],aspect_ratio=:equal);
Plots.scatter!([xs[d,i]],[ys[d,i]],[zs[d,i]], color = d,markersize = 6,label="Particula"*string(d),markershape=:star4);
else
Plots.plot!(xs[d,j:i],ys[d,j:i],zs[d,j:i],color=d,lw=2,label="",xlabel=L"X \mu m",ylabel=L"Y \mu m",bgcolor="black",alpha=a,xlim = [minimum(xs), maximum(xs)],ylim = [minimum(ys), maximum(ys)],aspect_ratio=:equal);
Plots.scatter!([xs[d,i]],[ys[d,i]],[zs[d,i]], color = d,markersize = 6,label="Particula"*string(d));
end
end
end;
gif(anim2d,"n-cuerpos.gif",fps=24)
end
Out[ ]:
animac (generic function with 1 method)
In [ ]:
function symN(N)
@variables X,Y,Z,M,x[1:N],y[1:N],z[1:N],v_x[1:N],v_y[1:N],v_z[1:N],t,m[1:N],a_x[1:N],a_y[1:N],a_z[1:N],r
r0=[X,Y,Z]
rs=[[x...]';[y...]';[z...]']'
vs=[[v_x...]';[v_y...]';[v_z...]']'
return x,y,z,v_x,v_y,v_z,r0,rs,vs,m,t
end
Out[ ]:
symN (generic function with 1 method)
In [ ]:
function ruido(S,I,R,p)
N=length(S)
ch=[1,0,-1]
for i in 1:N
r1=rand(ch,4).*p
S[i]=S[i]+r1[1]*(S[i])
#E[i]=E[i]+r1[2]*(E[i])
I[i]=I[i]+r1[3]*(I[i])
R[i]=R[i]+r1[4]*(R[i])
end
return S,I,R
end
ruido (generic function with 1 method)